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Abstract. - A non-perturbative analysis of the Bhatnagar-Gross-Krook (BGK) model kinetic 
equation for finite values of the Knudsen number is presented. This analysis indicates why 
discrete kinetic versions of the BGK equation, and notably the Lattice Boltzmann method, 
can provide semi-quantitative results also in the non-hydrodynamic, finite-Knudsen regime, up 
to Kn ~ 0(1). This may help the interpretation of recent Lattice Boltzmann simulations of 
microti ows, which show satisfactory agreement with continuum kinetic theory in the moderate- 
Knudsen regime. 



In the last decade, the lattice Boltzmann Method (LB) has developed into a very flexible 
and effective numerical method for the simulation of a large variety of complex flows, mostly in 
the macroscopic domain [1-3]. Fueled by relentless progress in micro, nano and bio-sciences, 
the recent years have witnessed a growing interest in exploring the possibility to enrich LB in 
the direction of describing micro-structured flows [4-6] . 

The LB method is based on a stylized stream-and-collide microscopic dynamics of fictitious 
particles, located on the nodes of discrete lattices and interacting according to local collision 
rules that drive the system towards a local equilibrium [2,3]. Mathematically: 

f z {x + c t At,t + At) - fi(2,t) = -LuAt [fifat) - ft 9 \x,t) 

where fi(x,t) is the probability to find a particle at position x and time t, moving along the 
lattice direction defined by the discrete speeds Cj (i = 1, b). The second term at the rhs of 
the above equations denotes relaxation towards a local equilibrium, the lattice analogue of a 
Maxwellian distribution in continuum kinetic theory: 

f^ 9 (x,t) = nwi{\ + (3ci -u + — \{ci ■ uf - u 2 ]. 

In the above, n(x, t) is the fluid density, u the flow speed, (3 = 1/c 2 is the inverse square speed 
of sound, and Wi a discrete set of weights normalized to unity. Finally, u, indicates a typical 
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time-scale relaxation (frequency relaxation) to local equilibrium, and governs the kinematic 
viscosity of the LB macroscopic fluid [2,7]. More general versions accounting for multiple-time 
relaxation [8-10] could be investigated but the simpler single-time (BGK) relaxation form will 
be sufficient for our present purposes. 

By taking suitable averages over molecular speeds, it can be shown that the resulting 
macroscopic quantities obey the Navier-Stokes equation of continuum mechanics. The fluid 
and current density are given by a linear superposition of the discrete distributions fi. 

b b 

n(x,t) = ^2fi(x,t), J(x,t) = ^2&ifi(x,t). (1) 

i=l i=l 

In order for those quantities to satisfy the exact hydrodynamic equations, it is required that 
the macroscopic fields do not vary appreciably on the scale of the mean free path A and the 
lattice spacing Ax. Based on the consolidated Chapman-Enskog background [11], one might 
be led to conclude that the range of applicability of LB methods is bounded by the domain of 
validity of the Chapman-Enskog method [12], the implication being that LB can only be used 
for strictly hydrodynamic purposes. Yet, such a restrictive stance is challenged by a number of 
recent numerical simulations [13-22] which clearly show that, by using appropriate boundary 
conditions, LB can reproduce some salient features of flows beyond the hydrodynamic regime, 
such as the onset of slip flow at finite-Knudsen numbers. 

In this Letter, we wish to propose a theoretical explanation for this rather unexpected valid- 
ity of LB in the beyond-Chapman-Enskog region. Our analysis is based on a non-perturbative 
solution of the lattice BGK (LBGK) equation, which overcomes the restrictions imposed to 
the standard Chapman-Enskog multiple scale analysis and extensions thereof [23,24]. 

Our analysis is confined to the bulk region of the flow, giving for granted that the use of 
proper boundary conditions is crucial to obtain correct results in actual LB simulations of 
microflows [25]. 

Let us refer for simplicity to the ld-continuum Boltzmann equation, written in BGK [26] 
form : 

d t f(x, v, t) + vd x f(x, v, t) = -w [f(x, v, t) - (x, t) 

where the local equilibrium f( eqS> (x,t) depends on x and t via the local velocity and density 
fields. Let us now consider the exact solution as provided, for each given velocity v, by an 
integration along the particle trajectory, from t to t + At: 



f(x + vAt, t + At)= e-" At f(x, t) + 
r At 

+u;e- uAt / e s "f( eg \x + vs,t + s)ds (2) 



Jo 

being At a generic time increment corresponding to the lattice time step. The above expression 
is exact, but purely formal, until one specifies a concrete procedure to compute the integral 
at the right hand side. This exhibits a quadratic functional dependence on f(x + vs, t + s), 
through the local equilibrium f eq . However, one can formally expand the integrand as: 

(x + vs,t + s) = J2 S -^f ieq) {x, t) (3) 

n=0 
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where D = dt + vd x denotes the streaming operator. By inserting this expansion into (J2J), we 
can formally solve the integral exactly, and obtain: 

f t+ . t = e~^f t+ e 1+ ° /u /<«> (4) 

where we have used the short-hand notation ft+ s = f(x + vs, t + s). This expression invites 
a number of useful considerations. Let us first recast it in the symbolic propagator form: 

ft+At = P(^At)f t + Q(ujAt, DAt)f ( t eq) (5) 

where 

P(cuAt) = e~ uAt (6) 
is the propagator from time t to time t + At, and 

p DAt _ p —uAt 

Q(u)At,DAt) = YT~Dju> (?) 

co-propagates the influence of the equilibrium at time t on the solution at time t + At. 
By neglecting the equilibrium variations on a scale At, i.e.: 

~ ti eq) 0<t'<At (8) 

the propagator e DAt acting on /( e 9) is such that 

e DAt f (e q ) = Q (Q) 

This leads to a much simplified version of the integral equation: 

ft+At = e-" M ft + (1 - e-" At )ti eq) (10) 

or, equivalently: 

ft+At - P("At)f t + Q LB (LuAt)ft q) (11) 
where we have defined the LB co-propagator as: 

Q LB (u;At) = (l-e-" At ). (12) 

The stick-to-equilibrium approximation @ is tantamount to retaining only the zero-th order 
term in the Taylor expansion @ , and consequently it is expected to work only for small values 
of At in the integral of (J2J ■ 

It is interesting to note that, already at this zero-th order level, and without any further 
assumption on ui, it is possible to identify a fully explicit discrete Boltzmann equation in BGK 
form. This reads: 

ft+At- ft = -UJLB{ft-f ( t eq) ) (13) 

with relaxation frequency lulb = (1 — e~ wA *). 

It is now instructive to analyze the two extreme limits uAt 3> 1 and LoAt «C 1, in the exact 
solution © first, and then in the discrete Lattice BGK solution (fill) . Let us begin with 
the former. Provided that the kinetic operator DAt remains bounded, we can identify the 
enslaving limit (uiAt 3> 1), characterized by: 

P(uAt) = e- uAt — » 
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DAt _ .-uAi 

Q(coAt,DAt) = - 6 >e DA *. 

1 + _L>/cj 

This is equivalent to enslave the solution to the local equilibrium, i.e. ft+At — ft+Af ^he 
opposite situation, (u>At <C 1), represents the free-molecular limit, in which one has: 

P(wAt) = e- uAt — ► 1 

DAt _ -uAi 

Q(uM,DM) = ' >0 

1 + D/uj 

corresponding to a free-flow (collisionless) solution ft+At — ft- 

A natural question arises on the nature of the same limits in the Lattice BGK equation 
The effect on the streaming propagator is the same, since the streaming term is inte- 
grated exactly in the lattice version. It remains to inspect the behaviour of Qlb{^>Al) in the 
two aforementioned limits. The enslaving limit (uAt 3> 1) yields: 

Q LB {ujAt) = l-e-" At — > 1 (wiB~0(l)) 

that is ft+At = ft^- I n view of the relation iJHJ, this is equivalent to state ft+At — ft+Af 
The opposite situation (u>At <C 1) yields: 

Qlb(wA£) = 1 — e~ wA * — ► {cj lb ^0). 

Thus, the two limits, full- enslaving and free-flow, are recovered by the LBGK, provided that 
the relaxation frequency is turned from bare lu to lulb- The question remains: what happens 
inbetween? A restrictive tenet is that Lattice BGK cannot work properly because the interme- 
diate regime involves all-order tensors, through the powers D n , which cannot be reproduced 
correctly in the discrete lattice because of lack of symmetry. On a more optimistic vein, one 
could counter-argue that since both extreme limits are correctly recovered, there might be 
hope that even inbetween the LB method could continue to provide a reasonable agreement 
with continuum kinetic theory. 

In order to analyze this point, we introduce the fine-scale Knudsen number Kn, defined 

as: 

Kn = X/S (14) 

where A = v/uj is the kinetic mean- free path and S is the smallest macroscopic length. 

Let us consider the most critical situation S = Ax — vAt, i.e. the macroscopic fields show 
appreciable variation on the scale of a single lattice site. Under these specific conditions, we 
obtain: 

uoAt = l/Kn (15) 

from which it follows that: 

P(uAt) = P{Kn) = e- x ' Kn (16) 

„DAt _ -1/Kn 

Q(wAi, DAt) = Q(Kn, DAt) = — — (17) 

1 + KnDAt 

Q LB {uAt, DAt) = Q LB {Kn, DAt) = 1 - e - 1/Kn . (18) 

Our Knudsen number-dependent approach is consistent with the fact that the full enslaving 
limit is recovered at Kn — * 0, which is equivalent to use a LBGK approach l|13|) with uj^gAt = 
(l — e~ uAt ) ~ 0{1). Starting from this Chapman-Enskog limit {Kn — > 0), it is then interesting 
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Fig. 1 - Plot of the percentage error E, as defined in 1201 . as a function of the Knudsen number up 
to the transition regime. Note that the error between the exact Boltzmann solution and the discrete 
form reaches a maximum of about ten percent at the top of the slip-flow regime (Kn ~ 0.1). 



to see whether the LB approach can be extended to the slip flow and transition regimes (up 
to Kn ~ O(l)) . In these regimes, the macroscopic fields fluctuate, in the most pessimistic 
case, on the same scale of the kinetic fields. Therefore, we may consistently assume that 
Df I f ~ v/S, and make the identification: 

DAt ~ Kn. (19) 

This is the zero-th order approximation relating DAt to Kn. More rigorously, one should 
consider a full series in Kn. However, to any practical purpose, this series can be truncated 
at the 0(Kn) without hampering the results up to Kn ~ 0(1). In this approximation, the 
propagators become: 

p Kn _ p -l/Kn 

Qlb(Ku) = 1 - e~'' Kn 

where Q(Kn) and QLs(Kn) represent the exact physical co-propagator effects and the lattice 
ones, in the range Kn < 1, respectively. The relative departure: 

E(Kn) = Qi^-Q^Kn) (2Q) 

is a quantitative measure of the spurious lattice-induced effects on the co-propagator. 

From figure l|T]l. where we report E{Kn) as a function of Kn, it is recognized that, up to 
Kn ~ 0.1, i.e. in the slip flow regime, the lattice Boltzmann approach does not differ from 
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the exact solution for more than to within a ten percent. Only as Kn exceeds 0.1, and the 
transition regime is entered, significant error build-up is observed. 

It is important to emphasize that the above figures apply to the worst-case scenario, in which 
macroscopic fields vary on the scale of a single lattice spacing (an extreme beyond-Chapman- 
Enskog situation). One may soften this assumption and assume that macroscopic fields exhibit 
significant variations on scales larger that a single grid spacing Ax, say: 

5 = hAx h = l,2... (21) 

h = 1 reproducing the previous worst-case scenario. This means that the Knudsen number 
is reduced by a factor h, Kn — > Kn/h, leading to a corresponding reduction in the error 
E(Kn). In addition, we observe that E(Kn) measures the error of the propagator, whereas 
physical observables result from the summation of discrete populations /, over the discrete 
speeds (see equations (|TJ) and consequently E(kn) as defined in l|20[l may well represent a 
pessimistic bound. In any event, the point of the present analysis is not to state the case 
for the accuracy of LB in the non-hydro dynamic regime, but only to point out that, even 
at finite-Knudsen, discreteness effects remain within fairly tolerable limits, at least for semi- 
quantitative purposes. Whether or not LB should be used instead of more accurate, and much 
more expensive, methods, such as Direct Simulation Monte Carlo, remains to be decided on 
a case-by-case basis. 

Summarizing, a non-perturbative analysis of the Boltzmann equation in BGK form, indi- 
cates that the LB method may continue to provide semi-quantitative agreement beyond the 
limits of the Chapman-Enskog theoretical framework. This could be of interest for the in- 
terpretation of Lattice Boltzmann simulations in the finite-Knudsen regime, including kinetic 
modeling of fluid turbulence [27] . Clearly, in order to realize the bulk properties highlighted 
by the present analysis, proper kinetic boundary conditions, [20,21], must be used in actual 
LB simulations of finite-Knudsen flows. 

R. Benzi, L. Biferale and F. Toschi are kindly acknowledged for useful comments and 
critical reading of this manuscript. 
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